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ABSTRACT 

The  computer  simulation  of  complex  systems  often  requires  efficient  recursive  algo¬ 
rithms  to  generate  various  noise  types.  Flicker  noise  often  occurs  in  frequency  and 
time  systems  and  long  samples  (a  million  or  more  lags)  are  needed.  In  June  of  1971 
two  papers  on  flicker  noise  simulation  were  published  independently.  Mandelbrot!1! 
based  the  simulation  on  the  sum  of  several  low  pass  filtered  white,  Gaussian  noises. 

He  showed  that  the  noise  level  and  the  filter  passband  of  each  noise  can  be  selected 
to  provide  an  approximate  flicker  noise  over  a  finite  but  arbitrary  spectral  range.  The 
precision  of  fit  to  the  flicker  spectrum  can  be  arbitrarily  good  depending  on  the  number 
of  independent  elements  used. 

Barnes  and  Jarvis!2l  based  their  simulation  on  a  cascade  of  lead-lag  filters.  Similarly, 
the  extent  and  goodness  of  fit  of  the  output  noise  depends  only  on  the  number  of  filter 
stages  used.  The  Barnes-Jarvis  procedure  has  an  exact  inverse  which  has  advantages 
over  the  Mandelbrot  method  in  some  statistical  applications  other  than  noise  simula¬ 
tion.  The  Barnes-Jarvis  method  can  be  expressed  as  an  ARIMA  model,  but  certain 
problems  arise  from  the  loss  of  significant  digits.  The  significant  digits  problem  can 
be  avoided  by  using  the  ARIMA  model  in  “factored  form”. 

Flicker  noise  is  a  noise  whose  power  spectral  density  (PSD)  varies  approximately  as 
the  reciprocal  of  the  Fourier  frequency  over  a  large,  but  finite  range.  Both  methods 
cited  can  be  used  to  simulate  other  power-law  noises  than  just  1  //  noise.  Any  noise 
with  a  PSD  that  varies  as  fa  for  a.  in  the  range  —  2  <  a  <  0  can  be  simulated  directly 
by  these  methods.  For  an  extension  of  the  range,  the  output  of  these  filters  can  be 
integrated  or  differentiated  to  obtain  any  power-law  desired.  The  Barnes-Jarvis  filter 
is  also  called  a  “constant  argument”  filter  since  the  output  phase  (45  deg.)  of  the  filter 
is  a  constant  relative  to  the  input  phase  for  a  sinusoidal  input  at  any  frequency.  A 
pure  integrator  has  a  90  degree  lag. 

THE  BARNES-JARVIS  FILTER 

Barnes  and  Jarvis  cascaded  a  series  of  lead-lag  filters  chosen  to  approximate  a  1//  spectrum. 
The  digital  equivalent  of  a  lead-lag  filter  can  be  expressed  with  a  recursive  filter!3!  in  the 
form: 


Yn  =  +  pn-  4>Pn-!  (i) 

where  Pn  is  the  input  to  the  filter  and  Yn  is  the  output.  The  coefficients,  4>  and  0,  determine 
the  poles  and  zeros  of  the  filter.  The  frequency  W  =  2irf  ( W  is  used  rather  than  u  to  aid  in 
comparison  to  computer  program  listings)  at  the  knee  in  the  transfer  function  (see  Fig.l) 
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determines  the  value  of  <f>  according  to  the  relation: 


<t>  =  l+  y(w0-^/wg  +  4)  (2) 

Conversely,  the  value  of  <f>  can  be  used  to  find  the  knee  frequency  according  to  the  inverse 
of  (2): 

w0  =  (1  -  *)/>/*  (3) 

The  same  functional  form  given  in  (2)  and  (3)  also  apply  to  the  9  coefficient  which  forms 
the  high  frequency  response  of  the  lead-lag  filter. 

One  next  selects  a  ratio,  R,  (R  >  l)  which  determines  the  frequency  ratios  of  successive 
filter  stages.  The  closer  that  R  is  to  unity  the  closer  the  filter  is  to  an  ideal  1/f  spectrum; 
but,  the  more  stages  are  required  to  fill  a  specific  frequency  range  (see  Fig. 2).  For  flicker 
noise  generation,  the  frequencies  corresponding  to  successive  knees  (determined  by  the  <£’s) 
grow  as  R2  as  do  the  frequencies  corresponding  to  the  9  coefficients.  In  each  stage,  the 
frequency  for  the  9  is  R  times  the  frequency  for  the  <f>. 

The  filter  stage  for  the  highest  frequency  filter  is  a  special  case.  First,  since  we  are  dealing 
with  real  functions  of  finite  length,  the  slope  of  the  spectrum  at  the  Nyquist  frequency 
(l/2T)  is  already  zero  of  necessity.  The  corresponding  9,  then,  can  be  taken  to  be  zero.  The 
next  problem  is  to  determine  a  frequency  corresponding  to  the  first  (highest  frequency)  <f>. 
This  selection  is  usually  done  by  trial  and  error  to  obtain  the  best  overall  approximation 
to  a  1/f  spectrum  (see  Fig.  2  and  Appendix  A).  The  first  <j>  is  typically  in  the  range  of  0.3 
to  0.5. 

A  simple  program  (in  BASIC)  to  determine  the  <£’s  and  0’s  follows: 


220  REM  COMPUTE  FACTORED  PHI’S  AND  THETA’S 
230  PH(1)=.45#:W=(1#-PH(1)) /SQR(PH(1)):TH(1)=0 

240  FOR  N=2  TO  M 
250  W=W/R 

260  TH(N)=l#+.5#*W*(W-SQR(W*W+4#)) 

270  W=W/R 

280  PH(N)=l#+.5#*W*(W-SQR(W*W+4#)) 

290  NEXT  N 


where  M  is  the  number  of  lead-lag  stages  to  be  used,  R  is  the  frequency  ratio,  W  is 
the  angular  frequencies  of  the  corresponding  ^’s  and  0’s,  and  .45  is  the  initial  <j>  selected 
from  trial  and  error.  The  symbol  indicates  double  precision  constants.  All  of  the 
calculations  should  be  carried  out  in  double  precision  except  for  storing  the  final  time 
series.  The  filter  equations  then  become: 


300  FOR  N=1  TO  NTOTAL 
350  Y(l)=PH(l)*Yl(l)+P 

360  FOR  1=2  TO  M 

370  Y(I)=PH(I)*Yl(I)+Y(I-l)-TH(I)*Yl(I-l) 

380  Y1(I-1)=Y(I-1) 

390  NEXT  I 

400  Y1(M)=Y(M) 

410  NEXT  N 


where  NTOTAL  is  the  total  number  of  points  to  be  generated.  The  inputs  to  the  filter 
are  the  P-values,  which  are  assumed  to  be  random,  normal  deviates  with  zero  mean  and 
a  variance  determined  by  the  particular  simulation  (unity,  here).  The  filter  output  is  the 
output  of  the  Mtfc,  or  final  stage:  Y(M).  Limited  storage  space  may  restrict  how  many 
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TRANSFER  FUNCTION 


R  =  S  FREQ.  RATIO  LOC(  FREQUENCY  ) 

N  =  4  FILTER  STAGES 

HK1)  =  .48  INITIAL  FHI  FIGURE  2a  FOUR-STAGE  BARNES/JARVIS  FILTER 


R  =  3.1  FREQ.  RATIO  LOG<  FREQUENCY  ) 

M  =  6  FILTER  STAGES 

PH(1)  =  .38  INITIAL  PHI  FI  (SIRE  2b  SIX-STAGE  BARNES/JARVIS  FILTER 


R  =  2.4  FREQ.  RATIO  LOG<  FREQUENCY  ) 

If  -  8  FILTER  STAGES 

PH(1>  =  .31  INITIAL  PHI  FIGURE  2c  EIGHT-STAGE  BARNES/ JARVIS  FILTER 


Figure  2.  Power  Spectral  Density  Times  Fourier  Frequency 
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values  can  be  stored,  but  the  filter  can  recursively  generate  as  many  values  as  desired  if 
all  the  values  do  not  need  to  be  retained  in  memory.  A  simple  and  reliable  algorithm  to 
convert  uniformly  distributed  random  numbers  on  (0,1)  to  normally  distributed  numbers 
with  zero  mean  and  unit  variance  is: 

310  P=0 

320  FOR  K=1  TO  6 

330  P=P+RND(1)-RND(1) 

340  NEXT  K 

where  each  pass  through  this  recursive  routine  produces  a  random,  normal  deviate. 

To  generate  a  million  values  would  overflow  the  storage  available  in  many  computers.  For 
the  simulations  reported  here,  every  100th  value  of  X(N)  (the  finite  integral  of  Y(M))  was 
stored  for  the  long  r-values.  For  the  smaller  r-values  only  parts  of  the  total  one  million 
data  points  were  used.  By  always  using  the  same  seed  for  the  random  noise  generator, 
one  can  extract  various  segments  of  the  same  noise  sample  without  having  all  data  points 
stored  at  one  time. 

The  filter  corresponding  to  lines  300  to  420  could,  in  principle,  be  replaced  with  an 
ARIMA(M,M-1)  model  by  eliminating  all  the  intermediate  F(i)  and  4T(r)  variables.  Un¬ 
fortunately,  such  a  procedure  requires  many  significant  digits,  and  even  double  precision 
is  totally  inadequate.  The  form  given  here  is  adequate  for  most  personal  computers,  but 
we  recommend  using  double  precision  anyway. 

Greenhalll4!  has  pointed  out  that  the  usual  convention  of  starting  out  a  filter  with  all 
past  values  set  to  zero  causes  systematic  errors  for  all  time.  In  effect,  the  flicker  filter 
has  significant  auto-correlations  for  all  lags  in  the  data  set  and  the  initial  zeros  produce 
a  non-representative  sample.  Greenhall  went  on  to  provide  an  initialization  algorithm 
which  circumvents  these  problems.  Appendix  B  gives  a  computer  program  to  calculate  the 
initialization  parameters. 

Figure  3  shows  the  Allan  variance  obtained  from  the  (finite)  integral  of  a  data  set  of  a 
million  flicker  simulated  numbers  using  the  algorithms  presented  in  this  paper.  That  is, 
we  simulated  the  phase  of  an  oscillator  perturbed  by  flicker  noise  modulating  the  oscillator 
frequency.  Since  the  available  computer  could  not  handle  a  million  points  all  at  once, 
sample  Allan  variances  were  obtained  by  storing  every  one-hundredth  phase  data  point, 
for  calculating  the  Allan  variances  for  the  larger  r-values.  For  flicker  FM,  the  Allan  variance 
is  constant!5!,  namely  /  x  S(f)  x  log(4),  where  S(f)  is  the  spectral  density  of  y  (Figs.  2  and 


THE  MANDELBROT  GENERATOR 

The  Mandelbrot  generator!1!  sums  a  set  of  independent,  band  limited  noises  to  approximate 
“Fractional  Brownian  Motions” .  This  can  be  realized  on  a  digital  computer  by  filtering  a 
white  noise  with  a  simple  low  pass  filter: 


Yn  =  4>oYn-i  +  Pn  (4) 

The  knee  in  the  Bode-plot  (Fig.  4)  for  the  A-th  low  pass  filter  is  located  at: 

("0'  (T^p  >  (5) 

That  is,  the  low  frequency  asymptote  of  the  filter  transfer  function  has  a  power  gain  of 
(1  -  4>o)~2  and  the  coefficient,  <£,  corresponding  to  the  frequency  at  the  knee  is  given  by 
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FIG.  4  BODE  PLOT  OF  LOW  PASS  FILTER 


Eq.(2),  above  (also  see  Appendix  A).  To  simulate  flicker  noise,  we  want  the  knees  of  the 
various  noises  to  fall  along  the  curve  5  =  1//. 

As  done  in  the  treatment  of  the  Barnes-Jarvis  model,  one  first  selects  a  frequency  ratio, 
R,  for  successive  noises,  and  then  calculates  the  corresponding  <f>  coefficient  using  Eq.  (2) 
for  each  noise  starting  from  the  first  frequency.  The  initial  (highest)  frequency  is  found  by 
trial  and  error  using  the  power  spectral  densities  as  calculated  in  Appendix  A.  Equating 
the  low-frequency  asymptote  of  each  noise  to  1//,  one  obtains  the  needed  variance  of  the 
input  white  noise  for  each  noise.  The  program  used  to  calculate  the  <t>’s  and  amplitudes, 
A(n),  is  shown  below: 


610 

620 

630 

640 

650 

660 

670 


FOR  N=1  TO  M 
W=C0/(RAN) 

PH(N)=l#+.5#*W*(W-SQR(W*W+4#)) 
A(N)=SQR(.5#*(l#-PH(N))*SQR(PH(N)J) 
P=0:FOR  J=1  TO  6:P=P+RND(l)-RND(l):NEXT  J 
Y(N)=P*A(N)/SQR(l#-PH(N)  A2) 

NEXT  N 


Lines  650  and  660  initialize  each  low  pass  filter  to  a  random  variable  whose  variance  is 
the  variance  of  its  steady  state  output.  This  takes  care  of  turn-on  transients  similar  to 
the  problems  in  the  Barnes-Jarvis  model.  The  coefficient  CO  is  chosen  by  observing  the 
theoretical  spectrum  at  the  higher  frequencies,  near  the  Nyquist  limit,  1/2T  (see  Fig.  5 
and  Appendix  A). 

The  program  which  generates  the  flicker  sample  is  shown  below: 

750  FOR  N=1  TO  NTOTAL 

760  Z=0 

770  FOR  1=1  TO  M 

780  P=0:FOR  J=1  TO  6:P=P+RND(1)-RND(1):NEXT  J 

790  Y(I)=Y(I)*PH(I)+A(I)*P 

800  Z=Z+Y(I) 

810  NEXT  I 

820  NEXT  N 

The  output  flicker  noise  is  Z  as  determined  recursively  for  each  value  of  N.  Figure  5  plots 
the  theoretical  power  spectral  density  times  /  on  a  log-log  plot  similar  to  Fig.  2. 


CONCLUSION 

Large  samples  (a  million  or  more  values)  of  pseudo-random  noise  which  approximates  a 
flicker  (or  1//)  noise  can  be  generated  by  rather  simple  recursive  functions.  The  range 
and  goodness  of  fit  can  be  selected  to  meet  any  specific  need  and  the  methods  are  not 
compromised  by  the  number  of  significant  digits  carried  on  most  computers.  The  Barnes- 
Jarvis  method  and  the  Mandelbrot  method  seem  to  provide  equally  good  results. 

In  regard  to  speed,  however,  the  Barnes-Jarvis  method  was  measured  to  be  over  4  times 
faster  than  the  Mandelbrot  method  running  six-stage  filters  of  comparable  performance. 
The  difference  in  speed  is  probably  due  to  the  fact  that  the  Mandelbrot  method  needs 
one  random  number  per  stage  per  point,  while  the  Barnes-Jarvis  method  requires  only 
one  random  number  per  point  regardless  of  the  number  of  stages.  The  Barnes-Jarvis 
method  ran  in  excess  of  300  points  per  second  using  a  compiled  BASIC  program  in  a 
micro-computer. 
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H  =  4  STAGES 
1HIT  UAL  :  18  C9 


FIGURE  5a  FOUR-STAGE  MANDELBROT  FILTER 


R  =  18  FREQ.  RATIO  LOG<  FREQUENCY  ) 

M  =  6  STAGES  „ 

f NIT  UAL  :  11  C8  FIGURE  5b  SIX-STAGE  MANDELBROT  FILTER 


R  =  5.8  FREQ.  RATIO  LOG(  FREQUENCY  > 

M  =  8  STAGES 

INIT  UAL  r  9  C8  FIGURE  5C  EIGHT-STAGE  MANDELBROT  FILTER 


Figure  5.  Power  Spectral  Density  Times  Fourier  Frequency 


Table  1.  ALF(ij) 


RATIO  =  2  PHI(l)  =  .3 
.31449 

.26035  .26333 
.11022  .28084  .27924 
.03183  .11429  .28500  .28480 
.00826  .03280  .11459  .28632  .28629 
.00208  .00850  .03278  .11468  .28668 
.00052  .00214  .00849  .03278  .11471 
.00013  .00054  .00214  .00848 .03278 
.00003  .00013  .00054  .00214  .00848 
.00001  .00003  .00013  .00054  .00214 

RATIO  =  2.5  PHI(l)  =  .325 
.34366 

.29249  .40792 
.07255  .29740  .44051 
.01254  .06936  .29817  .44641 
.00203  .01188  .06871  .29833  .44737 
.00033  .00192  .01175  .06861  .29836 
.00005  .00031  .00190  .01173  .06859 
.00001  .00005  .00030  .00190  .01173 
.00000  .00001  .00005  .00030  .00190 
.00000  .00000  .00001  .00005  .00030 

RATIO  =  3  PHI(l)  =  .35 
.37363 

.29450  .52478 
.04720  .28356  .55965 
.00548  .04290  .28232  .56381 
.00061  .00495  .04239  .28218  .56428 
.00007  .00055  .00489  .04233  .28217 
.00001  .00006  .00055  .00488  .04233 
.00000  .00001  .00006  .00055  .00488 
.00000  .00000  .00001  .00006  .00054 
.00000  .00000  .00000  .00001  .00006 

RATIO  =  3.5  PHI(l)  =  .375 
.40452 

.28370  .61773 
.03162  .26287  .64943 
.00265  .02792  .26111  .65214 
.00022  .00233  .02760  .26097  .65237 
.00002  .00019  .00231  .02757  .26096 
.00000  .00002  .00019  .00230  .02757 
.00000  .00000  .00002  .00019  .00230 
.00000  .00000  .00000  .00002  .00019 
.00000  .00000  .00000  .00000  .00002 


.28667 

.28676  .28676 
.11472  .28679  .28679 
.03278  .11472  .28679  .28679 
.00848  .03278  .11472  .28679  .28679 


.44752 

.29837  .44755 
.06859  .29837  .44755 
.01172  .06859  .29837  .44755 
.00190  .01172  .06859  .29837  .44755 


.56433 

.28217  .56433 
.04233  .28217  .56433 
.00488  .04233  .28217  .56433 
.00054  .00488  .04233  .28217  .56433 


.65238 

.26095  .65239 
.02757  .26095  .65239 
.00230  .02757  .26095  .65239 
.00019  .00230  .02757  .26095 .65239 
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RATIO  =  4  PHI(l)  =  .4 
.43644 

.26820  .69282 
.02195  .24240  .71998 
.00140  .01907  .24073  .72173 
.00009  .00121  .01888  .24062  .72184 
.00001  .00008  .00120  .01887  .24062 
.00000  .00000  .00008  .00120 .01887 
.00000  .00000  .00000  .00007  .00120 
.00000  .00000  .00000  .00000  .00007 
.00000  .00000  .00000  .00000  .00000 

RATIO  =  4.5  PHI(l)  =  .425 
.46951 

.25167  .75512 
.01575  .22407  .77784 
.00079  .01358  .22266  .77899 
.00004  .00068  .01347  .22259  .77905 
.00000  .00003  .00067  .01347  .22259 
.00000  .00000  .00003  .00067  .01347 
.00000  .00000  .00000  .00003  .00067 
.00000  .00000  .00000  .00000  .00003 
.00000  .00000  .00000  .00000  .00000 

RATIO  -  5  PHI(l)  =  .45 
.50390 

.23570  .80843 
.01164  .20820  .82727 
.00047  .01002  .20707  .82804 
.00002  .00040  .00996  .20702  .82807 
.00000  .00002  .00040  .00995  .20702 
.00000  .00000  .00002  .00040  .00995 
.00000  .00000  .00000  .00002  .00040 
.00000  .00000  .00000  .00000  .00002 
.00000  .00000  .00000  .00000  .00000 

RATIO  =  6  PHI(l)  =  .5 
.57735 

.20764  .89840 
.00685  .18306  .91130 
.00019  .00594  .18235  .91166 
.00001  .00017  .00591  .18234  .91167 
.00000  .00000  .00017  .00591  .18233 
.00000  .00000  .00000  .00017  .00591 
.00000  .00000  .00000  .00000  .00017 
.00000  .00000  .00000  .00000  .00000 
.00000  .00000  .00000  .00000  .00000 


.72185 

.24062  .72185 
.01887  .24062  .72185 
.00120  .01887  .24062  .72185 
.00007  .00120  .01887  .24062  .72185 


.77905 

.22259  .77905 
.01347  .22259  .77905 
.00067  .01347  .22259  .77905 
.00003  .00067  .01347  .22259 .77906 


.82807 

.20702  .82807 
.00995  .20702  .82807 
.00040  .00995  .20702  .82807 
.00002  .00040  .00995  .20702  .82809 


.91167 

.18233  .91167 
.00591  .18233  .91167 
.00017  .00591  .18234  .91168 
.00000  .00016  .00591  .18230  .91121 
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APPENDIX  A 


An  “exponential”  filter  or  simple  low  pass  filter  has  the  following  form: 


Yn  —  <f>Yn—i  +  an  (.41) 

where  0  <  <f>  <  1,  and  an  are  random  normal  deviates  with  zero  mean  and  variance  a2.  An 
equivalent  representation^  of  this  filter  in  terms  of  the  impulse  response  function  is: 

Yn  =  an  +  <f>an~i  +  ^2an_i  + . . .  (42) 

Taking  the  expectation  value  of  the  square  of  (A2)  one  obtains: 


ElYZ]  =  v2jr<f>2i, 


»= o 
i-<j>2 


(^3) 


since  the  an  are  independent.  The  initial  conditions  for  such  a  low  pass  filter  (e.g.,  as  used 
in  the  Mandelbrot  method),  can  be  taken  as  a  random  normal  deviate  with  this  variance. 
This  is  the  source  of  lines  650  and  660  in  the  text. 


Box  and  Jenkins  also  show  that  the  PSD  of  T(n)  is  given  by: 

glf\  —  _ _ 

1  +  <t>2  —  20cos(w) 


(A4) 


At  /  =  0,  the  value  S(0)  is  just  2a2/[(l  -  <j>0 )2]  and  we  define  the  cutoff  frequency  as  that 
frequency  for  which  S[f)  =  S(0)/2.  Using  the  first  two  terms  in  the  Taylor  series  expansion 
for  cos(2?r/)  in  (A4),  one  obtains: 

w0  =  (1  -  4>o)/V<h>  (-45) 


or,  equivalently:  _ 

<£o  =  1  +  (wo  -  \/w2  +  4) 
Also,  (see  lines  260,  280,  and  630): 

2cra  _  _L  —  V^O 

(1  —  <j)0)2  Wo  l  —  ^o 


(A6 


aa  =  2  V^o(!  -  to)  [A7) 

The  power  spectral  density  (PSD)  of  the  Barnes-Jarvis  filter  at  (angular)  frequency  W  can 
be  calculated  using  the  following  routine: 


500  S=2# 

510  FOR  J=1  TO  M 

520  S=S*(l#+TH(J)A2-2#*TH(J)*COS(W)) 

530  S=S/(l#+PH(J)A2-2#*PH(J)*COS(W)J 

540  NEXT  J 


and  for  the  Mandelbrot  method  the  PSD  can  be  calculated  using: 

800  S=0# 

810  FOR  J=1  TO  M 

820  S=S+2#*(A(J)A2)/(l#+PH(J)A2-2#*PH(J)*COS(W)) 

830  NEXT  J 

where  S  is  the  PSD.  S  is  multiplied  by  W/2k  for  Figures  2  and  5  to  display  the  errors  from 
a  perfect  1/f  noise.  The  PH(J)’s  and  TH(J)’s  are  those  calculated  for  the  appropriate 
noises. 
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APPENDIX  B 

Initialization  of  a  Barnes- Jarvis  Filter 


Following  GreenhallM ,  the  covariance  program  (in  BASIC)  can  be  written  in  the  form: 


1210 

1220 

1230 

1240 

1250 

1260 

1270 

1280 

1290 

1300 

1310 

1320 

1330 

1340 

1350 

1360 

1370 

1380 

1390 

1400 

1410 

1420 

1430 

1440 


’  COVARIANCE  PROGRAM 
FOR  N=1  TO  M 

C(1,N)=1#:D(1,N)=1# 

FOR  1=2  TO  M 

I1=I-1:F=ALF(I1)-BET(N) 

IF  IoN  THEN  F=F/(BET(I1)-BET(N)) 

C(I,N)=C(I1,N)*F 

F=ALF(Il)+BET(N)-ALF(Il)*BET(N) 
F=F/(BET(I1)+BET(N)-BET(I1)*BET(N)) 
D(I,N)=D(I1,N)/F 
NEXT  I 
NEXT  N 

FOR  1=1  TO  M 

FOR  J=1  TO  M 
SUM=0# 

FOR  N=1  TO  I 

F=C(I,N)*D(J,N)/(BET(J)+BET(N)-BET(J)*BET(N)) 
IF  IoN  THEN  F=F/(BET(I)-BET(N)) 

SUM=SUM+F 
NEXT  N 

RZ(I,J)=SUM*(ALF(I)-BET(I))*(ALF(J)-BET(J)) 

NEXT  J 
NEXT  I 


where  BET  (I)  =  1  -  PH  (I)  and  ALF(I)  =  1  -  TH(I).  The  desired  covariance  matrix  is 
R  Z(I,J)  • 


The  Choleski  square  root  of  the  covariance  matrix  can  be  computed  with  the  following 


BASIC  program: 

1450 

’  CHOLESKI  SQUARE  ROOT 

1460 

FOR  1=1  TO  M 

1470 

G=RZ(I,I) 

1480 

FOR  K=1  TO  1-1 

1490 

G=G-ALZ(I,K)*2 

1500 

NEXT  K 

1510 

ALZ(I,I)=SQR(G) 

1520 

FOR  J=I+1  TO  M 

1530 

TEMP=RZ(J,I) 

1540 

FOR  K=1  TO  1-1 

1550 

TEMP=TEMP-ALZ(: 

1560 

NEXT  K 

1570 

ALZ(J,I)=TEMP/ALZ(I,I) 

1560 

NEXT  J 

1570 

NEXT  I 

The  lower  triangular  matrix,  ALZ(J,I),  contains  the  coefficients  needed  to  initialize  the 
Yl(I)’s.  Table  B1  lists  the  coefficients  for  some  ratios  and  initial  PHI(l)  values.  For  fewer 
than  10  stages  the  matrix  can  simply  be  truncated. 


The  next  program  initializes  the  Barnes-Jarvis  filter: 
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1600  ’  INITIALIZATION 

1610  FOR  1=0  TO  M 

1620  P=0:FOR  J=1  TO  6:P=P+RND(l)-RND(l):NEXT  J 

1630  U(I)=P 

1640  NEXT  I 

1650  Y1(0)=U(0) 

1660 

1670  FOR  1=1  TO  M 

1680  Z(I)=0 

1690  FOR  J=1  TO  I 

1700  Z(I)=ALZ(I,J)*U(J)+Z(I) 

1710  NEXT  J 

1720  Y1(I)=Y1(I-1)+Z(I) 

1730  NEXT  I 

The  components  Yl(l),...,  Yl(M)  are  used  in  the  first  step  of  the  filter  (Lines  300  -  410). 
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